%% Learn P_q|r,s,t and p_q(dzeta_i), which gives state of each pixel (Lambda)
function [Lambda] = hmm(X)

% liczba iteracji
iterations = 5;

% wymiary obrazu
[M, N] = size(X);

% zainicjuj stany pikseli
Lambda = lambda(X);

% oblicz state transition probability
P = stp(Lambda);
% oblicz observation probability
p_q = op(Lambda, X);

% oblicz pseudo-likelihood
% L_prim = prod(prod(P)) .* prod(prod(p_q));
[tmp_c1, tmp_s1] = myprod(P);
[tmp_c2, tmp_s2] = myprod(p_q);

[Lc_prim, Ls_prim] = sci_multi(tmp_c1, tmp_s1, tmp_c2, tmp_s2);

% return;

%% iteruj do osiagniecia maksimum pseudo-likelihood 
for i = 1 : iterations
    % dla kazdego piksela...
    for m = 1 : M
        for n = 1 : N
            % zmien stan jednego piksela
            Lambda = ml(Lambda, m, n);
            Lc = Lc_prim;
            Ls = Ls_prim;

            % oblicz nowe state transition probability
            P = stp(Lambda);
            % oblicz nowe observation probability
            p_q = op(Lambda, X);
            % oblicz nowe pseudo-likelihood
            [tmp_c1, tmp_s1] = myprod(P);
            [tmp_c2, tmp_s2] = myprod(p_q);
            [Lc_prim, Ls_prim] = sci_multi(tmp_c1, tmp_s1, tmp_c2, tmp_s2);

            % przywroc stan piksela, jezeli nie ma poprawy
            [tmp_c, tmp_s] = mydiv(Lc_prim, Ls_prim, Lc, Ls);
            if comp(tmp_c, tmp_s, 1, 0)
                Lambda = ml(Lambda, m, n);
            end
        end
    end
end